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and {hjhjj^r)- These findings are compared to numerics and good agreement is found in the Mott 
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I. INTRODUCTION 

One of the most fascinating recent trends in physics 
of cold gases concerns atomic gases in optical lattices 
These systems offer "atomic Hubbard toolbox" 
Q that can be used for studies of condensed matter mod- 
els in a uniquely controlled manner. Cold atoms in op- 
tical lattices can be used for investigation of high Tc su- 
perconductivity 0, disordered systems 0, various spin 
models 0, novel quantum magnets P|, etc. Perhaps the 
most important version of the Hubbard model, that can 
be studied in optical lattices, is the Bose-Hubbard model 
(BHM) . This model is a prototypical system on which 
understanding of quantum phase transitions (QPTs) in 
boson systems is based 10], and it has important appli- 
cations in construction of a quantum computer 

Desp ite lots of theoretical studies on the BHM IT^ 
Q, E) E) 0, there is still lack of analytical predic- 
tions about some basic experimentally relevant quanti- 
ties. This motivates us to present here the results of a 
high order perturbative expansion in the tunneling cou- 
pling. These results, providing analytical approximations 
to different physical quantities in the Mott phase, should 
be helpful in the interpretation of experimental data. 

We focus on a one dimensional (ID) homogeneous sys- 
tem at average density of one atom per site, i.e., on the 
simplest BHM undergoing a QPT. We calculate the fol- 
lowing ground state quantities: energy per site £, atom- 
atom correlations C(r) = (a|aj+r), density-density cor- 
relations D{r) — (hjhj^r), and variance of on-site num- 
ber operator Var(n) = [(n^) — {nj)^]^/^. 

The Hamiltonian of interest, in terms of dimcnsionless 
variables used through the paper, reads 



M 



H = - J^(al+ia. + h.c.) + jYl " 

i=l i=l 

where we denote the number of atoms, and the number of 
lattice sites as M . Since physics of the BHM depends on 



M 



the J /U ratio only, we set [/ = 1 for convenience. With 
this choice, the critical point between the Mott insulator 
(MI) (0 < J < Jc) and superfluid (SF) (J > J^) is at 
Jc « 0.29 [3- 

Though the great deal of attention was recently de- 
voted to cold atoms in inhomogeneous lattices 0, IH 0| , 
the homogeneous systems described by the Hamiltonian 
can be realized in the nearest future in at least two 
setups. First, there is a recent experiment done in the 
Raizen group |l8j |. where a single, one dimensional, ho- 
mogeneous box is realized in a proper configuration of 
laser beams. After superposing a standing laser field on 
it, a ID homogeneous Bose-Hubbard model ^ with open 
boundary conditions can be achieved. Since it was al- 
ready demonstrated in that one can load this ID 
box with ultracold bosons and then count them very ef- 
ficiently, the studies of lattices with a desired number of 
atoms per site should be available (as discussed below 
on specific examples). Another experimental opportu- 
nity shows up after realization of the ring-shaped optical 
lattice proposed recently in Ref. 0. This time, a ID 
homogeneous lattice with periodic boundary conditions 
should be available for experimental investigations. 



II. THE METHOD 

Our findings come from a high order symbolic per- 
turbative expansion in the tunneling coupling. This 
method was first successfully applied to the calculation 
of a ground state and excited states of the BHM in one 
and two dimensional systems by Elstner and Monien . 
We compare perturbative expansions to numerical data 
obtained using the imaginary time evolution with the so- 
called Vidal's algorithm |23| , which is equivalent to Den- 
sity Matrix Renormalization Group scheme |2^. This 
allows for verification of accuracy of our analytical pre- 
dictions. The Vidal's algorithm calculations assume open 
boundary conditions, which breaks the translational in- 



2 



variance of the system. To minimize finite size effects 
during comparison between the numerics and perturba- 
tive expansions vaUd for infinite systems (where bound- 
ary conditions are irrelevant) we have calculated the cor- 
relation functions around the system center. 

We aim at calculation of high order perturbative cor- 
rections to different quantities of interest in the infinite 
Bose-Hubbard model. The perturbation theory is devel- 
oped around the Fock state |1, 1, . . .), where the numbers 
are boson on-site occupations. The expansion is done in 
the J parameter (^). 

In principle, the calculations can be performed by hand 
by perturbative determination of the wave- function up to 
a given order in the infinite system, and then subsequent 
calculation of expectation values in this wave-function. 
The attainable order of the expansion, however, is very 
limited (the wave-function can be determined up to the 
~ terms) so this method is not an option here. 

A better alternative is to perform a linked cluster ex- 
pansion (LCE) [221 that has been used so far in spin 
systems and the Bose-Hubbard model This 
method consists of two steps. First, one has to gener- 
ate all the clusters (sets of lattice sites in the BHM [Tsj l') 
that contribute to a given order of expansion. The largest 
of these clusters have the size comparable to the order of 
expansion. Then, one has to perform a perturbative ex- 
pansion in the relevant clusters and sum up the results 
properly. In the end, one gets perturbative expansion of 
different quantities valid for infinite system from analy- 
sis of finite clusters. Since the whole procedure can be 
implemented on a computer in a symbolic way, the high 
order expansions become feasible. 

The key for our calculations message from the LCE is 
the following: all the information about i-th order pertur- 
bative expansions in the infinite system is encoded in the 
small subsystems of the size ^ i. In accordance with this 
statement, we observe for every non-extensive observable 
that we consider, say O, that when we do a perturbative 
expansion 0\m = X^i'-^iU/^* [M is the system size) the 
following holds 



C{r) 



D{r) 



0^\M=0. 



M > M,(z), 



(2) 



i.e., the perturbative corrections become size indepen- 
dent for large enough systems. Naturally, Mdi) grows 
with i, but the growth is "reasonably" slow: Mc{i) ^ i 
(see Table IJ). We can not, of course, check expHcitly 
the relation ((SJ for arbitrarily large M, but it is clear 
from the LCE that the size independent expansion terms 
iOi\M>Ma) correspond to infinite system predictions. 

In our calculations we make a direct use of the relation 
(0) avoiding therefore generation of the cluster states. 
This simplifies the computer implementation of the whole 
procedure, but probably leads to more stringent require- 
ments on the computer resources. To be more specific, we 
fix the system size M, assume periodic boundary condi- 
tions to keep the system translationally invariant, and do 
a standard perturbative expansion leading to determina- 
tion of the ground state wave function up to a given order 



M >i + l M >i + l + r M >i + l and 2r < i 

TABLE L Conditions on the system size M for getting size 
independent perturbative predictions at i-th order to £, C{r) 
and D{r) - see These results are obtained for the system 
filled with one atom per site. 



(Appendix ^ . Having the ground state wave-function, 
we calculate perturbative corrections to the expectation 
values of different operators and study their dependence 
on the system size M. This way we easily get Mc{i) that 
guarantees lj2Jl: Tabled as well as the desired pertur- 
bative corrections valid for infinite systems. Notice that 
due to Q all finite size contributions are filtered out 
from the perturbative expansions for large enough M. 
To illustrate these findings we note that the 2nd order 
correction to the ground state energy per site in 3 sites 
and 3 atoms system is the same as in the infinite model 
at unit filling factor - a result that can be easily verified 
analytically. 



III. PERTURBATIVE EXPANSIONS 



To start, the ground state energy per site, S, satisfies 



58 _ 1257 44171 

9 81 1458 
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4902596 
6561 



-J 
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8020902135607 
2645395200 



J'^ + OiJ'*"). (3) 



Before proceeding further, it should be stressed that the 
fractions come from a symbolic calculation, whose details 
are presented in Appendix IXI 

Coming back to we notice that so far the largest 
published order of 8 expansion was the sixth fTpl . Nat- 
urally, the expansion of £ from jl5l] matches first three 
terms of Changing fractions into numbers one gets 
approximately the following sequence of non-zero coeffi- 
cients {-1, 1, 8,-16, 30, -747, -3032}, which shows that 
the series has a rather unpredictable form. 

To estimate accuracy of the expansion we plot (jSJ vs. 
numerical data for a fairly large system of 40 atoms in 
40 sites: Fig. ^ This plot shows that there is quite a 
good agreement between Q and numerics for J smaller 
than about 0.3, i.e., at least in a MI phase. The dis- 
crepancies present in Fig. ^ (and all other figures in this 
paper), may come from the following sources. First, our 
numerics is done in a finite system with open boundary 
conditions, while perturbative expansion yields results 
for infinite system. Second, we might need more expan- 
sion terms to get more accurate predictions. Third, the 
full perturbative expansion may fail to converge for large 
enough J's, probably J > Jc- In the worst case, the se- 
ries might be of asymptotic kind as was shown to be the 
case in some other systems i24j . Further discussion about 
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FIG. 1: The perturbative prediction for the ground state en- 
ergy per site from (|3J (sohd hne) vs. numerics for M — 40 
system (dots). 
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FIG. 2: The perturbative result for the variance of on-site 
number operator from Q (solid line) vs. numerics for A/ = 40 
system (dots). 



convergence of expansions presented here is beyond the 
scope of the present contribution. 

The variance of on-site number operator, Var(n), up 
to 0{J^^) terms satisfies 



Var(n) 



6979423019 
839808 



1441 
~36" 



32045 
648 
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3105413 



5184 



J' 
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5290790400 
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(4) 



which is illustrated in Fig. El As for the £ expansion, the 
agreement between perturbative prediction and numerics 
is good for J smaller than ~ 0.35. Thus, our result can 
be well applied to the system in the MI phase. Most 
interestingly, it shows that the variance of site occupation 
on the MI - SF boundary equals as much as 60% of its 
deep superfluid value in the limit of J — > cx) when the 
system is maximally delocalized (|B3ll . 

Expansion Q is useful, e.g., because there is an on- 
going experiment that aims at Var(n) determination in 
a ID untrapped setup |18| |. The measurement of Var(n) 
can be possible due to ability for a high efficiency single 
atom detection already shown in [l3|. It can be per- 
formed once extraction of atoms from a single lattice site 
will be demonstrated. Extracted atoms can be counted, 



and then all the remaining atoms from a lattice can be 
released and counted. Averaging the results of single 
site countings over the measurements where total num- 
ber of atoms is close to number of lattice sites, one should 
get experimentally Var(n). Another aspect of this exper- 
iment is that the lattice is blocked at ends with laser 
beams, which corresponds to open boundary conditions 
used in all our numerical calculations. 

Other quantities of interest are atom-atom correlation 
functions C(r) defined above. They were previously stud- 
ied numerically in the one-dimensional Bose-Hubbard 
model in p3| and perturbatively in two dimensional Bose- 
Hubbard models in |15| . They are important because an- 
other measurable quantity, the momentum distribution 
of atoms in a lattice, is expressed as ~ C'('') exp(iA:r) 
|25| , with k being atomic momentum. Here we list a few 
most important ones 



C(l) = 



C{2) 



C(3) = 



4 J - 8.r - 

39220768 , 
2187 



272 
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20272 
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-J' 



441710 
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94478400 

234862 



Jl3 + 0(jl5), 
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126040 - 
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Jl2 + 0( Jl4) 

7883333 
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220980576341 
1049760 ' 
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0(ji3) 



(5) 



The comparison between expansions ^ and numerics 
(Fig. 13) shows good agreement up to J equal to 0.25 ^ 
0.3, i.e., almost in an entire Ml phase. Additionally, © 
and our results for C(4 - • -7) indicate that C{r > 0) = 

o{.r). 

Expansions ((SJ reveal that atom-atom correlations 
take very substantial values at the critical point, e.g., 
C(l) at Jc equals about 0.8, i.e., 80% of its deep super- 
fluid value ljB3p : an interesting result showing that sys- 
tem wave function departs significantly from the |1, 1, . . .) 
state at the critical point. 

Finally we discuss the density-density correlations. 
They can be determined in a counting experiment almost 
the same as the one discussed for Var(n) measurement, 
except for the fact that now atoms from two depleted 
sites have to be counted first. A perturbation theory 
predicts l|njl with accuracy of 0{J^^). This data and our 
results for 0(4 ■■■7) reveal that D{r > 0) = 1 - 0{.P''). 

The comparison between the expansion and the nu- 
merics is presented in Fig. 0] there is a good agreement 
up to J equal to 0.2 ~ 0.225, which may be a little sur- 
prising concerning high order of One may attribute 
it to presence of open boundary conditions introducing 
small inhomogeneities in atoms density in our numerics. 
On the other hand, the effects of open boundaries will 
always be present in an experimental setup proposed in 
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FIG. 3: Atom-atom correlations as a function of the tun- 
neling J. Numerics for M = 40 system: C(l) (dots), C(2) 
(triangles), C(3) (x's). Lines: solid, dashed, dotted represent 
perturbative expansions of C(l), C(2), C(3), respectively. 

[Tsl. Thus, our numerical calculations might actually 
better represent the experimental situation than the ex- 
pansion ||HJ) derived for an idealized infinite lattice. 



Our perturbative predictions can be tested for self- 
consistency, e.g., the following identity can be derived 
from the eigenequation for the ground state energy: 

Using (j3Jl, Q and ||SJ one easily shows that indeed it 
exactly holds up to 0{J^^) terms, as expected from the 
orders of expansions presented in this paper. 

Though it is known that mean-field theory works badly 
in ID, it is instructive at this point to compare our find- 
ings to predictions of the most popular mean- field, i.e., a 
Gutzwillcr variational wave function P, H^ . This ap- 
proach predicts in ID, at unit filling factor, that the 
ground state is |1, 1, . . .) for < J < 1/11.6, which is in- 
terpreted as a Ml phase. It implies that in this J range: 
£■ = 0, Var(n) = 0, C{r) = 0, D{r) = 1. Our results 
analytically quantify the amount of discrepancy between 
these mean- field predictions and the exact ones. It would 
be also instructive to compare our findings to predictions 
of the perturbatively improved Gutzwiller approach . 



AT^ 1^ r4 76 150638 8 4897282 415922848153 ,3 1022120948444278027 ,4 

D(l) - 1-4J +—J + 14696640 ^ 7777461888000 ' 

_ 100 4 2128 6 1156462 8 6848011 ^ 10808763042127 5150051155340205251 ,4 

^ ' ^ ^ ~ ^ ~3~ 243 729 ^ 44089920 3888730944000 ' 

13064 3727066 . 1588041877 ^ 1710030328933 2208787916976404357 ,4 
n ']) = 1 7° J 7° J 71^ J 714 

^ ' 27 243 8748 2755620 370355328000 



(6) 



IV. SYMMETRY OF PERTURBATIVE 
EXPANSIONS 

Expansions of £, C(r), and D{r) have well defined par- 
ity with respect to the J — J transformation. That 
is explained below, and provides an insight into the 
J ^ — J symmetry of the BHM. Additionally, it helps 
in determining the accuracy of our perturbative predic- 
tions. Indeed, the symmetry of perturbative expansions 
proven below implies that if the last calculated nonva- 
nishing term shows up in the r-th order, the expansion 
term in the r + I order vanishes. 

First, the ground state energy per site, Eq. satis- 
fies £{J) = £{—J). To see it one performs the canonical 
transformation 

a, ^ {~iyh, (7) 

which is a bosonic equivalent of the Shiba transformation 
used to prove a similar symmetry property of the Fermi- 
Hubbard model: Sec. 2.2.4 of 27]. After the transfor- 



mation, Hamiltonian ^ reads 

M M 

H ^ J ^(a^^^ifli + h.c.) + jYl ~ 

iff (i) the system is infinite; (ii) the system consists of 
even number of sites and periodic boundary conditions 
are applied; (iii) open boundary conditions are chosen. 
To the end of this paper we assume that one of these 
conditions holds, and that there exists a unique ground 
state. The transformed Hamiltonian is a J — > — J ver- 
sion of (Q), so we get the desired prediction about the 
symmetry of £{J)- 

Second, we focus on C(r)|j>o = (+|a]aj-|_r|+), where 

l+> = E^{".}(«Ir---(«l)""|vac) (8) 
{»<} 

is a ground state of ^ with J > [C^nj are some 
coefficients and dijvac) = 0]. Using {Tj), one gets that 
the ground state of the Hamiltonian with J < is 
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FIG. 4: Density-density correlations as a function of the tun- 
neling J. Numerics for M = 40 system: Dil) (dots), D{2) 
(triangles), -D(3) (x's). Perturbative expansions @: solid, 
dashed and dotted lines represent -D(l), D{2), D{3), respec- 
tively. 

(up to a phase factor): 

i-> = E(-i)^'"^^{".}(«Ir---(4/r"ivac), 

where C^m} are the same as in (jSJ. An exphcit 
calculation shows that C(r)|j<o = (— |ajaj-|_r|— ) = 
{~iyC{r)\j>o, therefore C(r) ^ {-iyC{r) as J ^ - J: 
see (5) for an example when ^ = M oo. 

Third, following above calculation, one shows that 
D{r) and Var(n) does not change under J sign inversion. 
The latter result is shadowed by explicit assumption of 
J > during Taylor expansion leading to Q). 

Finally, we note that the ground state of two compo- 
nent Bose-Einstein condensate in an optical lattice was 
found to change after inversion of J sign 28] . Our paper 
provides tools for a closer inspection of this interesting 
finding. 



V. SUMMARY 

In this paper, we have derived and discussed analytical 
predictions on important quantities describing physics of 
the Mott phase: variance of on-site number operator, 
atom-atom and density-density correlation functions. We 
have also explored, for the first time to our knowledge, 
the fundamental J ^ — J symmetry of the BHM. We 
expect that our findings can be useful in interpretation of 
ongoing experiments performed in homogeneous lattices 
filled with ultracold atoms. 
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J.Z. was supported by Polish Government scientific funds 
(2005-2008). 



APPENDIX A: PERTURBATIVE EXPANSION 

Here we sketch some of technical details of the nonde- 
generate perturbative expansions that we perform. 

To introduce the notation: the Hamiltonian (^) is re- 
written as Hq + XV, where A = — J, Hq = ^ Si^i 

l),l> = E.^i(alia^+h.c.). 

The ground state wave-function and eigenenergy up to 
a given order are 

5 S+1 

|v&)=EAn^(r)), E=J2x^£r, (Al) 
lvM>=E«r'l*)' Ho\^)=e,\^), m=5,,. (A2) 

i 

To keep the wave function ljAl|l normalized to unity we 
must satisfy = 1 + 0{X^+^). It implies that 

{ip{0)\ip{Q))' ^ 1 and m 

1 '■"^ 

(^(O)l^(r)) =--E(^(»)|^(r-»)), r^l,...,S. (A3) 

1=1 

The wave function HA1(I follows the eigen-equation 

iHo + XV)\^) ^ E\^), (A4) 

which in the A° order implies that |<^(0)) = ao"''|0) = |0) 
(after proper choice of basis), and £q ~ cq. In higher 

orders Aq'^' = (0|(p(r)) is found from ljA3|l . 

The coefficients a^^^^ l|A2p are obtained from projec- 
tion of HA4|I onto basis states {|j 7^ 0)} HA2p 

«r(e, - eo) + E(^l^b>r'^ = gf.af-'^). (A5) 
j k=i 

The eigenenergies are found from a projection of (|A4|I 
onto |0) [23 

j fe=i 

where r = 1, . . . , 5 -I- 1. Equations HA3|) . (|X5l) and (|X6|) 

have to be solved order by order starting from r = 1. 
In our calculations 

1^(0)) :=|0) = |1,1,...), (A7) 

and 

\i) ^\nu,n2i,...,nMi), ^ ^'^{uki - l)nki. (A8) 

Efficient manipulations on Fock states (|A8|I are cru- 
cial for getting high order expansions on a com- 
puter. We have achieved it by using a fast mapping 
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{nii, n2i, . . . , n Mi } i which was provided by the hash 
functions - see [Sy for their general description. 

Since ljA8|l is an integer number and the action 
of V onto a Fock state produces factors of the form 
^ynk±i + ly^^ife (jik is the atom occupation of the k-th 
lattice site), we see from (|A5|) and (IA6|) that every a'f-' 
coefRcient will be in general equal to 



bi{r,i},s) 
^ c({r,i},s)^ 



E 



l[pi{r,i},s,k), (A9) 



k=l 



where b{{r,i},s) and c{{r,i},s) are integer num- 
bers and p{{r,i},s,k) are prime numbers satisfying 
p{{r, i}, s,k) < r + 1 - notice that when we develop the 
expansion around the Fock state ljA7p the highest on-site 
occupation in the r-th order is r -I- 1 because it originates 
from r consecutive actions of V operator onto (|A7(I . Sim- 
ilarly, one finds that kaia.x{{r,i}, s) < r. 

We have done the symbolic calculations by writing a 
proper C code. Each set of {p{{r, i}, s,k) : k = 1 . . . fcmax} 
has been encoded in bits of a single integer number, while 
every b{{r,i}, s)/c{{r,i}, s) ratio has been allocated as 
a Gnu Multi Precision rational number !3ll|. The first 
choice minimizes the memory requirements on storage 
of the prime number part, while the latter one allows 
for avoiding overflows occurring at high order expansions 
when b{{r,i},s) and c({r, j},s) are allocated as integer 
numbers. To reduce usage of memory the Smax param- 
eter has to be set as small as possible during allocation 
process. By running the program we found that Smax = 2 
is sufficient for all the calculations performed in this pa- 
per. 

The results of the symbolic calculations were directly 
verified by calculation of the perturbative expansion in a 
numeric (as opposed to a symbolic) way. There we treat 

(r) 

al as double precision numbers avoiding complications 
of symbolic manipulations. Results obtained numerically 
are in agreement with the symbolic ones within accuracy 
provided by the double precision format. 



APPENDIX B: NONINTERACTING SYSTEM 

It is instructive to briefly summarize here the predic- 
tions for various observables in the limit of negligible 



atom interactions, say U = case. 

The N- particle ground state of the non-interacting 
M- site Hubbard model with periodic boundary con- 
ditions is given by 



\SF) = 



1 



M 



N 



vac) 



(Bl) 



where aj|vac) — and {SF\SF) — I - notice that (|Bip is 
nothing else than the N- particle zero momentum ground 
state • The ground state eigenenergy equals —2JN. 
After some algebra (|Blll can be rewritten to the form 



\SF) 



m 



E 



|ni, . 



, nM) 



(B2) 



where the sum goes over all sets of boson on-site oc- 
cupation (integer) numbers {rii} such that X^i"* = 
N and < < TV. The wave-function ||B2|) 

allows for an easy calculation of different expecta- 
tion values. Indeed, the summation can be ex- 
pressed as J2{n,} = J2n^=0 ' ' ' E^m=0 Sn^ + -+nMM, thcn 

the Kronecker delta can be represented as i5/^„i = 
^ J^^ d(p[exp{iip)]'-~"^ , and finally the integration over 
if can be replaced by the integration over z — exp(ic/?) 
on the unit circle surrounding the origin of the complex 
z plane. Usage of the residue theorem allows for getting 
that 



Var(n).V--^, = W - 



for any integer N > and M < r. 

In this paper we are interested in the case of N/M — 1 
and M > 1, i.e.. 



Var(n) = 1, C(r) = 1, D{r) = 1. 



(B3) 
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